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Abstract 

| /' | Quadrature formulas for spheres, the rotation group, and other compact, homogeneous man- 

ifolds are important in a number of applications and have been the subject of recent research. 
The main purpose of this paper is to study coordinate independent quadrature (or cubature) 
formulas associated with certain classes of positive definite and conditionally positive definite 
kernels that are invariant under the group action of the homogeneous manifold. In particular, 
we show that these formulas are accurate - optimally so in many cases -, and stable under 
an increasing number of nodes and in the presence of noise, provided the set X of quadrature 



> 
in 



nodes is quasi-uniform. The stability results are new in all cases. In addition, we may use 
these quadrature formulas to obtain similar formulas for manifolds diffeomorphic to S™, oblate 
spheroids for instance. The weights are obtained by solving a single linear system. For § 2 , and 
the restricted thin plate spline kernel r 2 logr, these weights can be computed for two-thirds of 
a million nodes, using a preconditioned iterative technique introduced by us. 

1 Introduction 

Quadrature formulas for spheres, the rotation group, and other compact, homogeneous manifolds 
are important in many applications and have been the subject of recent research [l6][l7][2§[3l][32]^ 



The main purpose of this paper is to study quadrature (or cubature) formulas associated with 
certain classes of positive definite and conditionally positive definite kernels that are invariant 
under the group action of the homogeneous manifold, and to show that these formulas are accurate 
and stable, provided the set X of quadrature nodes is quasi-uniform. The invariance of the kernels 
^ is a key ingredient in giving simple, easy to construct linear systems that determine the weights. 

The weights themselves have size comparable to 1/N, where N is the number nodes in X. For S 2 , 
and the restricted thin plate spline kernel r 2 logr, these weights can be computed for very large 
values of N using the preconditioned iterative technique described in [IT] . 

The quadrature formulas developed here are for the general setting of a compact, homogeneous, 
n-dimensional manifold M that is equipped with a group invariant Riemannian metric g%j and its 
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associated invariant measure dfj,(x) = \J det(gij(x))dx. From the point of view of current applica- 
tions, the two most important of these manifolds are § 2 and SO(3). However, other homogeneous 
spaces, such as Stiefel and Grassmann manifolds, are arising in applications fll[7j , and we expect 
that in the future our results will be applied to them. The type of quadrature formula we will be 
concerned with here has the form 



/ f{x)dn{x) = V <*/(£) =: Q(f), f G C(M), (1.1) 
JM 77Z 



where the set X C M of centers/nodes is finite. The weights {c^}^x are chosen so that quadrature 
operator Q integrates exactly a given finite dimensional space of continuous functions, V. 



In the case of S n and SO (3), popular choices for V are spaces of spherical harmonics 16,26 



31 32 , and Wigner D-functions [17 . Work also has been done on compact two-point homogeneous 



manifolds |6| and on general compact homogeneous manifolds 34 , with V being chosen to be a 
set of eigenfunctions of a manifold's Laplace-Beltrami operator. We will use the term polynomial 
quadrature for methods that integrate such spaces exactly; this is because, for spheres and a few 
other spaces, V consists of restrictions of harmonic polynomials. 

The quadrature methods developed here are kernel methods; they use a space V consisting 
of linear combinations of kernels. In [40] , Sommariva and Womersley used spaces of rotationally 
invariant radial basis functions (RBFs) and spherical basis functions (SBFs) to derive a linear 
system of equations for the weights eg. However, neither accuracy, control over the size of the 



weights, nor stability was addressed in 40 . In this paper these and other issues are dealt with 
employing recent results developed by us in [TTJ[l9j[20] . 

Accuracy is measured in terms of the mesh norm h, which is defined in section [2| All the kernels 
we deal with are associated with a Sobolev space W™, for some m in N. Previous error estimates 



for quadrature with positive definite kernels on § 2 were given in 26 . For the general case of a 
homogeneous manifold, the error we obtained is 0(h m ), for functions in W™. However, on S n or 
SO(3), we get even better rates, in fact optimal: if a function is in C 2m , then the order is 0(h 2m ). 
For example, the thin-plate spline kernel restricted to S 2 has m = 2, so, for a function in C 4 (S 2 ), 
the error would be 0(h ). 

In the course of studying accuracy for Q in various spaces, we obtain new error estimates for 
interpolation/reproduction problems on two-point homogeneous manifolds. First, for a general such 
manifold, the estimates hold for / G W^"(M), which is the "native space" of the kernel. Second, 
for spheres (new when reproduction is required) and real projective spaces, the estimates apply to 
/ G WJf , with re/2 < \i < m. Thus the estimates allow an "escape" from the native space, in the 
sense that they are for functions not smooth enough to be in that space. 

These quadrature formulas are stable both under an increase in the number of points and in the 
presence of noise. If the number of points is increased, then the norm of the quadrature operator 
remains uniformly bounded, as long as the level of quasi-uniformity is maintained. Thus there is 
no oscillatory "Runge phenomenon." To examine the effect of noise, we assume the measured func- 
tion values differ from the actual ones by independent, identically distributed, zero mean random 
variables. Under these conditions, the standard deviation of the quadrature formula decreases as 

To illustrate how the method works, consider a positive definite SBF kernel <f>(x ■ y), x,y G 
S n , and let V = span{0(x • £)}gex- We want to integrate s G V exactly; that is, we require 
X^ex c £ s (£) = J§™ s(x)dfi(x). Doing so results in a systems of linear equations for the weights. If 
A is the interpolation matrix with entries = ■ ry), £, 77 G X, then the vector of weights c 
satisfies (Ac)^ = L n 4>(£, ■ x)dfi(x). At this point, we encounter an apparent difficulty. We have to 
compute every one of the integrals L n 4>{£ • x)dfx(x) for every £ £ X. 
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The rotational invariance of the SBF allows us to overcome this difficulty. Because of rotational 
invariance, all of the integrals are independent of £, and thus have the same value Jo; the system 
then becomes (Ac)^ = Jo- The constant Jo only needs to be computed once for a given kernel. 
(For many SBFs/RBFs, values of Jo are known; see [40] .) The same is true for group invariant 
kernels on M, as we will see in section 2.2 below. The point to be emphasized here is that group 
invariance of the kernels allows us to deal with quadrature as an interpolation problem. Without 
it, the problem requires computing many integrals and in fact becomes prohibitively expensive, 
computationally. 

Numerical tests of the quadrature formulas were carried out for M = § 2 , in connection with 
the SBF kernel $(x • y) = (1 — x ■ y) log(l — x ■ y), x, y € §> 2 . This kernel is the thin-plate spline 
r 2 log r restricted to the sphere, and it corresponds to one of the Sobolev spaces mentioned above, 
namely, W 2 (§ 2 )- For these tests, the sets of quasi-uniform nodes were generated via three different, 
commonly used methods: icosahedral, Fibonacci (or phyllotaxis), and quasi-minimum energy. The 
number of nodes employed varied over a substantial range, from a few thousand to two-thirds of 
a million. Weights corresponding to these nodes were computed using a pre-conditioning method 



developed in 11 . The tests themselves focused on the accuracy and stability of the method, both 
in terms of increasing the number of nodes and adding noise. The tests, which are discussed in 
section [5j gave excellent results, in agreement with the theory. 

There are situations where the manifold M involved is not a homogeneous space, but quadrature 
formulas can still be obtained. If M is diffeomorphic to a homogeneous space, then it is possible 
to obtain quadrature weights for M from the ones for the corresponding homogeneous space. In 
section [6j we will show how this can be done for E> n . We will then apply this to the specific 
case where M is an oblate spheroid (e.g., earth with flattening accounted for), which is of course 
diffeomorphic to S 2 . 

The paper is organized this way. Section [2] begins with a brief discussion of positive defi- 



nite/conditionally positive definite kernels, notation, and, in section 2.1 interpolation. Section 2.2 



contains a derivation and discussion of kernel quadrature formulas, with special emphasis on the 
role played by group invariance of the kernel employed in the formula. In section [2~3} the questions 
of accuracy and stability mentioned in the introduction are taken up. The results obtained there 
are aimed at invariant kernels, such as Sobolev and polyharmonic kernels discussed in sections [3] 
and|U 

Sobolev kernels on a compact manifold are positive definite reproducing kernels for the Sobolev 
space W™, m > n/2. In section[3j we study these in terms of their invariance, interpolation errors, 
and properties of their Lagrange functions and Lebesgue constants. Finally, in section 3.4 we look 
at their use in quadrature formulas. 

Section [4] is devoted to a very important class of kernels on a compact, two-point homoge- 
neous manifold: the polyharmonic kernels. These kernels, which may be either positive definite 
or conditionally positive definite, are Green's functions for operators that are polynomials in the 
the Laplace-Beltrami operator. On spheres, they include restricted surface splines, and on SO(3) 
similar kernels. All of these are given in terms of simple, explicit formulas. The whole section is a 
self-contained discussion of these kernels, culminating in their application to quadrature formulas. 

The results from various numerical tests that we conducted are discussed in detail in section [5j 
Finally, in section [6] we discuss ways of using quadrature weights for a compact, homogeneous 
manifold M to obtain invariant, coordinate independent weights for manifolds diffeomorphic to M. 
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2 Interpolation and Quadrature via Kernels 



The spaces that we will work with will be homogeneous manifolds, ultimately. However, for inter- 
polation we need very little in the way of structure. In fact, we could take our underlying space M 
to be a metric space. 

The set X C M will be assumed finite. Its mesh norm (or fill distance) h := sup^gM dist(x, X) 
measures the density of X in M, while the separation radius q := | inf^ e x dist(£, £) determines 

the spacing of X. The mesh ratio p := h/q measures the uniformity of the distribution of X in M. 

We say that a continuous kernel k : M x M — > R is (strictly) positive definite on M if, for 
every finite subset X C M, the matrix A with entries A^^ := K^,rj), £,77 G X, is positive definite. 
Conditionally (strictly) positive definite kernels are defined with respect to a finite dimensional 
space II := spanj^fc : M — > M.}™ =1 , where the t/Vs are linearly independent, continuous functions 
on M. In addition, given a finite set of centers X C M, where we let N := j^X be the cardinality 
of X, we say that II is unisolvent on X if the only function ip G II for which ip\x = is y = 0. 
This means that {tpk\x}kLi is a linearly independent set in K . . Given that II is unisolvent on X, 
we say that the kernel k is conditionally positive definite if for every nonzero set {a^ G Mj^gx such 
that ^£ o^V'fc(C) = 0, k = 1, . . . ,m, one has 

a^a v n(C,rj) > 0. (2.1) 



2.1 Interpolation 

Positive definite and conditionally positive definite kernels can be used to interpolate a continuous 
function / : M — > M, given the data f\xi by means of a function of the form 

m 

s = o^«(-, C) + ^ frfcV'fc, where = 0, k = 1, . . . ,m. (2.2) 

£gx fc=i ^eJf 

We will denote the space of such functions by Vx . 

In the case where the kernels are RBFs or SBFs, the space II is usually taken to be either the 
polynomials or spherical harmonics with degree less than some fixed number. 

We now turn to the interpolation problem. Let = ipk\x, k = 1, . . . , m, and define the N x m 
matrix ^ = Vt^ • • • ^ m ]- In addition, let a = (a^)^ e x and b = (pi ■■■b Tn ) T . The constraint 
condition that ^tfl^fcCO = can now be stated as ^ T a = 0. Requiring that s interpolate a 
function / G C(M) on X is then f\x = s\x = Aa + \ti>, where Ag„ = k(£,T]). Written in matrix 
form, the interpolation equations are 



A * W a\_f f\x j (2 3) 

Jmxl 



* T Or. 



Using the constraint condition ^> a = and the positivity condition (2.1), one can easily show 
that the matrix A on the left above is invertible. In addition, the interpolation process reproduces 
II; that is, if / G II, then s = f. Finally, much is known about how well s fits /. In many cases, 
the approximation is excellent (cf. [45] and references therein). 

In the sequel, we will need the Lagrange function centered at £ G X, x$ G Vx- We define X( to 
be the unique function in Vx that satisfies x^iv) = that is, X£ ls 1 when x = £ and when 
x = 7] G X, 7] / £. 
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2.2 Quadrature 

We will now develop our quadrature formula for a C°°, n- dimensional Riemannian manifold M that 
is a homogeneous space for a Lie group Q 144] . This just means that Q acts transitively on M: for 
two points x, y G M there is a 7 G 5 such that y = 73;. Equivalently, M is a left coset of Q for a 
closed subgroup. 

S 2 is a homogeneous space for SO (3). In fact, the Lie group Q is a homogeneous space for itself. 
Thus, SO (3) is its own homogeneous space. Homogeneous spaces also include Stiefel manifolds, 
Grassmann manifolds an many others. All spheres and projective spaces belong to a special class 
known as two-point homogeneous spaces. Such spaces are characterized by the property that if two 
pairs of points x,y and x',y' satisfy dist(a;,y) = dist(x',y / ), then there is a group element 7 G Q 
such that x' = 72; and y' = *yy. While spheres and projective spaces are compact, there are also 
non compact spaces: W 1 and certain hyperbolic spaces also belong to the class. 

Invariance under Q plays an important here. We will assume throughout that the kernel n is 
invariant |43[ §1.3.4] under Q. Moreover, we will take d/x to be the (/-invariant measure associated 
with the Riemannian metric tensor for M [43[ §1.2.3]. Thus, for all x,y G M, 7 6 Q, and / 6 Li(M) 
we have 

K(-yx,-yy) = K(x,y) and / f(x)dfi(x) = f (jx)dn(x) . (2.4) 



In particular, all SBF kernels on S n , which have the form <p(x ■ y), are invariant, as is the standard 
measure on S n . The following lemma is a consequence of k and dfi being invariant. 

Lemma 2.1. The integral J(y) '■= f M K(x,y)dfi(x) is independent ofy. 

Proof. Take z € M to be fixed. Because of the group action on M, we may find 7 G <5 such that z = 
■jy. Prom the invariance of n under 7, we have n(x, y) = k(jx, z). Hence, J(y) = J M K("fx, z)d[i{x). 
However, the integral being invariant under Q then yields 

J(y) = / ^(72;, z)d^(x) = I k(x, z)d^(x) = J(z), 
Jm Jm 

which completes the proof. □ 

Since J(y) is independent of y, we may drop y and denote it by Jo, which we will do throughout 
the sequel. Note that Jo may be 0. In addition, we define these quantities. 

J k : = fu^k(x)dfi(x), k = l,...,m; 



J :— (Jl • • • J n 

1 := 



We point out that, for the constant function 1 on M, 1 = l\x is the column vector in ~R N with all 
entries equal to 1. 

The result below gives a formula for the integral of a function s G Vx and provides a system of 
equations that determine the quadrature weights. It follows the one given in [401, with rotational 



invariance replaced by the invariance under Q proved in Lemma 2.1 



Proposition 2.2. Let k and dfi satisfy \2.J$ and suppose that c and d are the N x 1 and m x 1 

column vectors that uniquely solve the (N + m) x (N + m) system of equations 

Ac + ^d = J 1 and f T c = J. (2.5) 
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If s G Vx j then 



Jm 

In addition, if \i £ Vx ^ e Lagrange function centered at £, i/ien we a/so /iaue 

JM 

Proof. Integrating s(x) from ( |2.2| ) results in this chain of equations: 

/si 7T1 n 

s(x)dfi(x) = y^o^ / K(x,£,)dn(x) + y^&fc / ^k{x)dii{x) 
a 7~y Jm r~. Jm 

£,eX s / fc=l s / 



(2.6) 



(2.7) 



(2. 



J l T a + J T b 



(Jol T J T ) 



Jo 

a 
b 



Using (2.3), with f\x replaced by s\x, together with the invertibility and self adjointness of A, we 
obtain 

a \ = ( A' 1 ( Jq1 \\ T ( s|x 
b 1 { J { mx i 



T T T\ 



T i 

c s\x- 



Multiplying 



by A and writing out the equations for c, d yields the system (2.5). Combining 



the two previous equations then yields (2.6). Moreover, from (2.6), with s replaced by X£ and the 
values of six replaced by those of X£ on the set X, we obtain the formula for in (2.7). Finally, 
the bound on the right in (2.7) follows immediately from the integral formula for eg. □ 



As we mentioned earlier, the quadrature formula for / G C(M) is obtained by replacing / with 
its interpolant in Vx- To that end, we define the linear functional Qv x (f) that will play the role 
of our quadrature operator. 

Definition 2.3. Let f G C(M) and let sf G Vx be the unique interpolant for f , so that Sf\x = fix- 
Then, we define the linear functional Qv x '■ C(M) — > M. via 



Qv x (f)-= s f (x)dfi{x) = V c ? /(£), 
Jm fex 



(2.9) 



where the c^'s, the weights, are given in Proposition 2.2 

The invariance assumption on the kernel k produces a system with the same attractive feature 
as one for an SBF <f>(x ■ y). The integral Jq depends only on k and M; it is entirely independent of 



X. This is also true of the other Jfc's. As a result the equations (2.5) defining the weight vector 



c only depend on X through function evaluations. The integrals Jo, J\, ■ ■ ■ , J m are all known in 
advance and are independent of X. 

It is important to note that, without the invariance of k, obtaining the system for c would require 
computing integrals of the form J* M k(x , £)d[i(x) for each £ G X. This follows just by looking at 
equation (2.8) in the derivation of (2.5). This would make finding c numerically very expensive, 



not only because each integral would have to be computed for every £ G X, but also because the 
whole set would have to be recomputed whenever X was changed. 
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2.2.1 Weights 



In many cases, the weights appearing in the quadrature formula above may be interpreted as coming 
from simple interpolation problems. In particular, if k is strictly positive definite and II = {0}, 
then the system (2.5) becomes Ac = Jq\. This is the set of equations for interpolating the function 
f(x) = 1/Jq. To obtain the interpolation problem for quadrature with k being merely conditionally 
positive definite with respect to II, start with ty T c = J. This equation completely determines the 
orthogonal projection of c, cm := Pc, onto the range of ^: The standard normal equations give 
C || = Pc = ^(^ T ^) _1 ^ T c = ^(xj/^vl/)- 1 J. Thus, C|| is known. Next, let c± := P ± c, which 
obviously satisfies Pc± = 0, so ty T c± = 0, and, consequently, the following system: 

Ac ± + ^d = J 1 - A^i^^y 1 J and V T c± = 0. (2.10) 

This is the interpolation problem (2.3), with fx = </ol — A 1 ^( 1 ^ Tl ^)^ 1 J. The final weights are then 

c = c± + y(^ T ^)- 1 J. (2.11) 



Note that (2.10) may be solved for cj_, without also having to solve for d as follows. Start by 
eliminating d from (2.10). Multiply both sides of (2.10) by P L . Since P ± c± = c± and P- 1 ^ = 0, 
we get 

P ± AP ± c± = JqP^I - P ± A^(y T yy 1 J. (2.12) 

Because k is conditionally positive definite relative to II, P AP is positive definite on the or- 
thogonal complement of the range of ^. Restricted to that space, it is invertible. Carrying out the 
inverse gives us c±. 

Often, the space II contains the constant function; that is, 1 € II. When this happens, the term 



with Jo drops out of (2.12), which then becomes 

P ± AP ± c± = -P ± Ac {l = -P ± A^(^ T ^y 1 J. 

This equation is homogeneous in A and is therefore independent of Jo- This implies that c± can 
be determined independently of Jo- 

There is another consequence of having 1 G II. First, the column vector 1 is in the range of ^f, 



so PI = 1. Let (cj_) be the average of cj_. Since Pc± = 0, we have N(c±) 



l T c ± 



(PlYc ± = 0. 



Second, since 1 S II, we also have that the quadrature formula is exact for it, and so Q(l) 
vol(M) = l T c = N{c). Using c = c± +cm, along with a little algebra, yields 



<ci>=0, (c) = ( C ||) 



vol(M) 



provided 1 € II. 



(2.13) 



Signs of weights in quadrature formulas are usually nonnegative. This is true for the polyno- 
mial quadrature formulas developed for spheres [3l] and for two-point homogeneous manifolds (6j. 
Numerical experiments by Sommariva and Womersley |40| produced negative weights for kernels 
formed by restricting Gaussians to S 2 . On the other hand, their experiments involving the thin-plate 
splines (cf. section [4]) restricted to S 2 resulted in only positive weights. 

Determining for what kernels and with what restrictions on X all of the weights are positive is 
an open problem. 
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2.3 Accuracy and Stability of Q 



V x 



There are two important questions that arise concerning the quadrature method that we have been 
discussing. First, how accurate is it? Second, how stable is it? 

The accuracy of Qy x depends on how well the underlying space of functions, Vx m our case, 
reproduces functions from the class to be integrated. Specifically, we have the following standard 
estimate: 

\Qv x (f) - [ f(x)dfi(x)\ < \\s f - f\\ Ll{u) < (voKM)) 1 / 2 !!^ - /|| t2(M) . (2.14) 
Jm 

This inequality converts estimates of the accuracy of Qv x (/) to error bounds for kernel interpola- 
tion. These are known in many cases of importance. 

Stability is related to how well the quadrature formula performs under the presence of noise. 
Lack of stability can amplify the effect that noise in f\x will have on the value of Qv x (/)• Stability 
also relates to performance as the number of data sites increases. Standard one-dimensional equally 
spaced quadrature formulas that reproduce polynomials can be quite unstable, due to the well- 
known Runge phenomenon. 

A measure of stability is the C(M) norm of the quadrature operator. From the definition of 



Qv x i n (2.9) and equation (2.7), we have that 



\\Qv x \\c{M) = M - Yl llxdUi(M)- ( 2 - 15 ) 

When M is a compact manifold, this bound can be given in terms of the Lebesgue constant, 
Ay x = maXa-gMiX^ex ^he reason is that 

\\Qv x \\cW) ^ E HxdUi(M) < J2 I \x^)W{x) < vol(M)Ay x . (2.16) 

Let us briefly see how noise affects Qv x (/)• Suppose that at each of the sites £ we measure /(£) + 
f£, where is a zero mean random variable. Further, we suppose that the v^s are independent 
and identically distributed, with variance a 2 ,. Our quadrature formula gives us Qv x (f + u ) rather 
than Qv x (f)- Because the v^s have zero mean, we have that the mean E{Qy x (/ + v)} = Qv x (f)- 
The variance of Qy x (f + v) is thus 

°Q = E{{Q Vx (f + u)- Q Vx (f)) 2 } = E{Q Vx (u) 2 }. (2.17) 
We can both calculate and estimate ctq: 

Proposition 2.4. Let be independent, identically distributed, zero mean random variables having 
standard deviation o~ v . Then the standard deviation aQ satisfies 

°Q = °l\\ c \\l < o-^llcllillclloo < vol(M)cj 2 Ay x max||^||L 1 (M) (2-18) 



Proof. We begin by evaluating the term on the right in (2.17). to do this, we need to compute 



E^^Ujj}. Since the u^s are i.i.d., we have E{u^u v } = a 2 5^^. It follows that o~q = E{Qy x (v) 2 } 



o"^||c||2- Moreover, ||c||2 < ||c||oo||c||i. Combining this with (2.7) and the previous equation results 



in (2.18). □ 
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3 Sobolev Kernels 



There are two classes of kernels that we will discuss here. The first class was introduced in 19, § 3.3], 
in the context of a compact n-dimensional Riemannian manifolcQM equipped with a metric g. This 
class comprises positive definite reproducing kernels for the Sobolev spaces W^iM), as defined 
in [3 24 . For M a homogeneous space with g being the invariant metric that M inherits from a Lie 



group Q, we will show below that these kernels are invariant under the action of Q. 



The second class, which was introduced and studied in 20 , comprises polyharmonic kernels on 
two-point manifolds - spheres, projective spaces, which include SO(3), along with a few others. 
These kernels are conditionally positive definite with respect to finite dimensional subspaces of 
eigenfunctions of the Laplace-Beltrami operator and they are invariant under appropriate transfor- 
mations. We will discuss these in section [4] below. 

The Sobolev space W^M), to G N, is defined as follows. Let (•, -) g>x be the inner product for 
a Riemannian metric g defined on TM. X , the tangent space at x S M. This inner product can also 
be applied to spaces of tensors at x. We denote by V fc the k th order covariant derivative associated 
with the metric g, and let d\i be the measure associated with g. For W^M), define the inner 
product 

m „ 

(/, h) mM : = (/, h) wm(M) : = V / (V k f, V k h) ax d/i(x), (3.1) 

k=o Jm 

and norm ||/||^ M := (/, /)m,Mi where /, h : M — > R are assumed smooth enough for their W™ 
norms to be finite. The advantage of this definition is that it yields Sobolev spaces that are co- 
ordinate independent and can also be defined on measurable regions Vt C M. Using the Sobolev 
embedding theorem for manifolds |3, §2.7], one can show that if to > n/2 these spaces are repro- 
ducing kernel Hilbert spaces, with K m being the unique, strictly positive definite reproducing kernel 
for W 2 m (M); that is, 

f( x ) = (/(■) 5 «m( a: >-))m,M 

In the remainder of this section we will discuss invariance, interpolation error estimates, Lagrange 
functions, Lebesgue constants, and quadrature formulas derived from K m . 

3.1 Invariance of n m 

We now turn to a discussion of the invariance of n m under the action of a diffeomorphism that is 
also an isometry. We will then apply this to the case of a homogeneous space. Here is what we will 
need. 

Proposition 3.1. Let M be a compact Riemannian manifold of dimension n with metric g. If <3? : 
M — > M is a diffeomorphism that is also an isometry, then the kernel K m satisfies K m ($(x), 3>(y)) = 
K m (x, y) and dfj,(<&(x)) = dfi(x). 

Proof. The proof proceeds in two steps. The first is showing this. Let / : M — > R and let /* = /o$ 
be the pullback of / by Then, 

(V k f*,V k h*) gx = (V k f,V k h) gMx y (3.2) 



We will follow a technique used in [25| Proposition 2.4, p. 246] and in 20 . Let (it, 4>) be a local chart, 
with coordinates v? ' = 4> 3 (x), j = 1, . . . , n for x £ 51. Since $ is a diffeomorphism, ($(iX), <f> o <3> -1 ) 



l See 

etc. 



19 



for a discussion of the Riemannian geometry involved, including metrics, tensors, covariant derivatives, 
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is also a local chart. Let ifi = (f>o <I> _1 , and use the coordinates = ip^{y) for y € ^(il)- The choice 
of coordinates has the effect of assigning the same point in W 1 to x and y, provided y = &(x) - 
i.e., v?(x) = v 3 {y). Thus, relative to these coordinates the map $ is the identity, and consequently, 
the two tangent vectors (gfj)j/ € T y M. and {~^j)x £ T X M. are related via 

_d_\ =d$ 

So far, we have only used the fact that $ is a diffeomorphism. The map $ being in addition an 
isometry then implies that 

— — \ = ld<$> (—) d$ ( — )) = / — — \ 
dvi , dv k / Hx) \ x \^v?) , X \du k )/ Hx) \du^ du k / x ' 

The expressions on the left and right are the metric tensors at y = <J>(x) and x; the equation implies 
that, as functions of v and u, gjk(v) = gjk(u). From this it follows that the expressions for the 
Christoffel symbols, covariant derivatives and various expressions formed from them also will be 
the same, as functions of local coordinates. In addition, note that the local forms for /* and 
at u = u(x) are /* o cj)~ l and o 4>~ 1 , respectively, and those for / and h at v = v(y) are / o if;" 1 
and h a tp" 1 . At u = u(x) and v = v(y), y = $(a;), / o ij)~ l (v) = / o $ o = /* o <f)~ l (v), 

and similarly for h. Again, this is functional equality in local coordinates, so matching partial 



derivatives are also equal. Consequently, (3.2) holds. 

The equality of the coordinate forms of the metric gjk, which was established above, implies 
invariance of the Riemannian measure: 

dfi(x) = ^det(g jk (u))d n u = ^det(g jk (v))d n v = dfj,($(x)). (3.3) 

Using this we see that 

/ (V fc /*,V fc ^> dn(x) = f (V k f,V k h) d^x) 
Jm Jm 

= [ (V k f,V k h) d^(x)) 
Jm y y ' 



f (V k f,V k h) duty) 
Jm y ' y 



Finally, from this and (3.1), we have invariance of the Sobolev inner product: 

(/ $ ,^ $ )m,M = (f,h)m t M (3.4) 

The second step is to show that the kernel is invariant. Since K m is a reproducing kernel, we 
have f(x) = (/(■), Km(x, -))m,M- By (3.4), we also have f(x) = (/*(■), K m (x, $(-)))m,M- Replacing 



x by <J>(x) then yields 

/(*(s)) = /*(*) = (/*(•), Km(*{x), H-))) mM . 

Finally, replacing / by /* 1 above gives us 

f(x) = (/(•), K m ($(z),$(-))>m,M, 

from which it follows that K m (&(x), &(y)) is also a reproducing kernel for W™(M). But, reproducing 
kernels are unique, so K m ($(x), = n m (x,y). Thus K m is invariant. □ 
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Homogeneous spaces have two properties that allow us to use the proposition just proved. First, 
they inherit a Riemannian metric invariant under the action of the Lie group Q. Second, the action 



of a group element produces an isometric diffeomorphism 25 43 . These observations then yield 
this: 

Corollary 3.2. Let M be a homogeneous space for a Lie group Q and suppose that M is equipped 
with the invariant metric g from Q . Then the reproducing kernel K m for W™(M.) is invariant under 
the action of 7 £ Q. 



3.2 Error estimates for interpolation via K m 

Recall that the accuracy of the quadrature formula associated with n m is directly dependent on 
error estimates for interpolation via n m . To obtain these, we will first state a theorem that provides 
estimates on functions with many zeros, quasi-uniformly distributed over a compact manifold. 



Theorem 3.3 ( |20[ Corollary A. 13] "Zeros Lemma"). Suppose that M is a C°° , compact, n- 
dimensional manifold, that 1 < p < 00, m £ N, and also that u £ W™(M). Assume that m > n/p 
when p > 1, and m >n when p = 1. Then there are constants Co = Cq(M.) and C\ = C\{m, k,M.) 
such that if u\x = and X C M has mesh norm h < Co/m 2 , then 

IMIw*(M) - C\h m fc ||u||wm( M ). (3.5) 

Using this "zeros lemma" we are able obtain an estimate for \\sf — /||l 2 (m)> provided / £ W™(M) 
and sj is the interpolant for /. 

Proposition 3.4. Let m > n/2, f € W™(M), and let sj be the K m -interpolant for f from Vx- 
Then, with the notation from Theorem 3.3, if h < Co/m 2 , we have 



s f ~ /IIl 2 (m) < Ci^ m ||/||w 2 m (M)- (3-6) 



Proof. Clearly Sf — f G W™ and (sf — f)\x = f\x — fx = 0. Applying Theorem 3.3 to Sf — / with 
k = yields \\sf — /||l 2 (m) < Cih m \\sf — f\\w^-(M)- Since the space W^M) is the reproducing 
Hilbert space for n m , the interpolant Sf minimizes \\g — /||vy 2 m (M) among all g £ W^M). Thus, 



taking g = yields \\s f - f\\ W ™{M) < II - /||w 2 m (M) = ll/llvF 2 m (M), from which the inequality ( [3^6 ) 
follows immediately. □ 



The bounds in (3.6) hold whether or not M is homogeneous. In one respect, the bounds are not 
as strong as we would like. The assumption that / is in the reproducing kernel space W™ precludes 
estimates for less smooth /, in W2 (M), with k < m. Stronger bounds do hold in important cases, 
though. For example, in section [4J we will see that the stronger result for / G W?? holds for the class 
of "polyharmonic" kernels, where M can be a sphere or some other two-point homogeneous space. 
We conjecture that, at least for K m , the stronger result holds for general C°° compact, Riemannian 
manifolds. 



3.3 Lagrange functions and Lebesgue constants 

The Lagrange functions {x^}^ex associated with K m have remarkable properties. For quasi-uniform 
sets of centers, Lebesgue constants are bounded independently of the number of points and the L\ 
norms of the Lagrange functions are nicely controlled. Specifically, we have the result below, which 
holds for a general C°° metric g. 
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Proposition 3.5 ( |19[ Theorem 4.6] |18[ Proposition 3.6]). Let M be a compact Riemannian 
manifold of dimension n, and assume m > n/2. For a quasi-uniform set X C M, with mesh ratio 
h/q < p, there exist constants Cm and Cm,™ such that if h < Cm/™ 2 , then the Lebesgue constant 
Ay x = max-rgM X^ex IXfO^)! associated with n m satisfies 

Av x < CM,mP n - 

In addition, we have this uniform bound on the L\-norm of the Lagrange functions: 



max||x € || Ll ( B 



[) < Cp jmi ] 



Proof. In the statement of the theorem used above, we have made explicit the p dependence of the 



bound on Ay x found in the proof of 19 Theorem 4.6]. We have also adapted the notation there to 
that used here. The bound on ||xcllii(M) follows from |18j Proposition 3.6], with s = x§ and p = 1. 



In that proposition, the coefficients A P:71 
and so ||xelUi(M) 



q n / p 5p, 71 and so \\A P 



n/p 



For p = 1, H^li, 



^ Cp y mQ n - 



□ 



3.4 Kernel quadrature via k Ti 



The positive definite Sobolev kernel K m is invariant under the action of £/, by Corollary 3.2 This 



is the bare minimum requirement for a kernel to be able to give rise to a computable quadrature 
formula. The other properties of n m established in the previous sections give us the result below 
concerning accuracy and stability. 

Theorem 3.6. Let X CM. be a finite set having mesh norm px < P- Take A = (ft m (£, r]))\^^x and 
suppose f £ WrpCM). Then the vector of weights in Qy x is c = J$A~ l l, the error for Qy x satisfies 
\Qv x (f) ~ Jm/^H - ^ m H/llw 2 ™(M); and the norm of Qy x is bounded by ||Qy x ||c(M) < C M ,raP n ■ 
Finally, the standard deviation defined in Proposition 2.4 satisfies the bound oq < C ' p, m ,MO~uh n ^ 2 ■ 



Proof. The formula for the weights is a consequence of two things. First, by Corollary |3 . 2| the kernel 
is invariant and so by Lemma 2.1 the integral J*^jK m (x, ■)dp{x) is a constant, namely Jq. Second, 



since the kernel is positive definite, Proposition |2.2| provides the desired formula for the weights. 



The error estimate is a consequence of (2.14) and Proposition 3.6, and the norm estimate follows 



from (2.16) and Proposition 



3.5 



Finally, the bound on <tq is a consequence of Proposition 
Proposition 3.5, and of the fact that p n q n = h n . 



2.4 



□ 



There are two significant implications of this result. The first is just what was noted in see- 
namely, the weights are obtained by directly solving a linear system of equations. The 



tion 2.2 



second is that the measures of accuracy and stability hold for any X with mesh ratio less than a 
fixed p. 

There are several drawbacks. In the case of a general homogeneous space M, the formulas for 
kernels n m are not yet explicitly known. This may be less of a problem when specific cases come up 
in applications. Also, the error estimates, which provide the chief measure of accuracy, hold for / 
in the space ^^(M). We would like to "escape" from this native space (reproducing kernel space) 
and have estimates for less smooth /. As we shall see in the next section, the situation is much 
improved in the case of spheres, projective spaces, and other two-point homogeneous manifolds. In 
the case of S 2 , not only do we have the requisite kernels, but in important cases we can also give 
fast algorithms for obtaining the weights (cf. section [5]). 
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4 Polyharmonic Kernels 



For spheres, 50(3), and other two-point homogeneous spaces (cf. [25j pgs. 167 & 177] for a list), 
one can use polyharmonic kernels, which are related to Green's functions for certain differential 



operators. These include restrictions of thin-plate splines 20 , which are useful because they are 
given via explicit formulas. Many of these kernels are conditionally positive definite, rather than 
positive definite. 

The differential operators are polynomials in the Laplace-Beltrami operator. Since, on a com- 
pact Riemannian manifold —A is a self adjoint operator with a countable sequence of nonnegative 
eigenvalues Xj < Xj+i having +oo as the only accumulation point, we can express a polyhar- 
monic kernel in terms of the associated eigenfunctions — A0 JiS = Xj(j)j )S} s = 1, . . . ,dj, dj being the 
multiplicity of Xj. To make this clear, we will need some notation. 

Let m £ N such that m > n/2 and let ir m £ Il m (IR) be of the form ir m (x) = J2T=o c ^ xl ' i where 
c m > and let C m be the 2m-order differential operator given by C m = vr m (— A). We define JcN 
to be a finite set that includes all j for which the eigenvalue ir m {Xj) of C m satisfies n m (Xj) < 0. (In 
addition to this finite set, J may also include a finite number j's for which 7r m (Xj) > 0.) We say 
that the kernel k : M x M — > R is polyharmonic if it has the eigenfunction expansion 

k(z> y) ■= Yl R i ( Yl faAtffcAv)) ' where k i '■= \ * m . ^ ' J - f J J i 4 - 1 ) 

j =0 V s=1 / [ arbitrary, j t J . 

The kernel k is conditionally positive definite with respect to the space Hj = span{0, s : j £ 
3 'j 1 < s < (ij}. In the sequel, we will need eigenspaces, Laplace-Beltrami operators and so on 
for compact two-point homogenous manifolds. Without further comment, we will use results taken 
from the excellent summary given in [6]. This paper gives original references for these results. 

Polyharmonic kernels are special cases of zonal functions, which are invariant kernels that 
depend only on the geodesic distance d(x, y) between the variables x and y; d is normalized so 
that the diameter of M is ir — i.e., < d(x,y) < ir. To see that polyharmonic kernels are zonal 
functions, we need the addition theorem on two-point homogeneous manifolds. If P^ a '^ denotes 
the I th de gree Jacobi polynomial, normalized so that P^ a '^\l) = r(€+i")r(a+i) > then this theorem 
states that for any orthonormal basis {cf>j tS }^ J =1 of the eigenspace corresponding to Xj we have 

( u f \ P^^l ( \W + /3)r(/3 + i)r(£+f +/3) 
> (bj s {x)(pn s(v) = c^r^A cosle dix.y))), cp := ; -rr — ; , 

~[ m, yy> 3 £] y v v ,ynh r(/3 + ^)r(£ + /? + i) 

where the parameters e and f3, and of course Xj, depend on the manifold. These are listed in the 
table below. In summary, every polyharmonic kernel has the form n(x,y) = Q(cos(e~ 1 d(x,y))), 
where 

00 -2 

$(t) = Y ZjCejP*^®, -1<*<1- (4.2) 

3=0 

An important class of polyharmonic kernels on S n are the surface (thin-plate) splines restricted 
to the sphere. The surface splines on W a are defined in |45[ Section 8.3]; their computed 
in [11 . Both kernel and Kj are given below in terms of the parameter s = m + n/2; they are also 
normalized to conventions used here (see also [30]). The first formula is used for odd n and the 
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Manifold 


e 


P 


•\j 


S n 


1 


n—2 
2 


i(i + n - 1) 


MP n 


2 


n-2 
2 


2j( 2 j + n-l) 


CP n 


1 





i(i + f) 


Quaternion P™ 


1 


1 


i(i + ^) 


Caley P ie 


1 


3 


JXJ + 12) 



Table 1: Parameters for compact, two-point homogeneous manifolds . 



second for even n. These kernels are conditionally positive definite of order m. 

( (-l) s +i(l - £) s , s = m _| G N+i 

(-iy+\l-tyiog(l-t), s = m-%eN. \ ( 1.3) 

Kj = Cs,n r Q ( +7+ n) , where, for even n, j > s. 
where the factor C Sin is given by 

c s , n : = r^T( s + i)r (s + - +2 ' +2 

2 [ 1, sGN. 

For the group SO(3) (= MP 3 ), the following family of polyharmonic kernels, which are similar to 
those above, was discussed in [21] : 

k(x, y)=( sin ( ^L^l) ) , m >2. (4.4) 



Here uj(z) is the rotational angle of z E SO(3). Adjusting the normalization given in |2Tj Lemma 
2] to that used here, the polynomial ir(x) associated with k is tt(x) = Y\™=oi x + 1 — 4i/ 2 ). We 
remark that derivation of the kernels given in 21 used the theory of group representations along 
with a generalized addition formula that holds for any compact homogeneous manifold (see [12]). 

Remark 4.1. For two-point homogeneous manifolds, the Sobolev kernels K m are in fact polyhar- 
monic. The reason is that K m is the Green's function for the linear operator C m = Ylk=o J ^ k 
However, by [2o| Lemmas 4.2 & 4.3], C m = 7r m (— A), where vr m (x) = x m + Yl™=o c v x>/ ' ■> so ^ 
has to have the form (4.1). Also, it is easy to show that 7r m (Aj) > for j = 0, 1, . . ., and that 
K m ,i = vr m (Aj) _1 ~ Xj m for j large. 

There are two reproducing kernel Hilbert spaces associated with a polyharmonic kernel k. We 
will take Kj > for j £ J . These spaces are often called "native spaces," and are denoted by M K 
and J\f K , j ■ Consider the inner products 

</,,). = y, E ^ ^ (f-^j = E E ^ 

j=0 s =l 3 jf£J s=l 3 

The first of these is the inner product of the reproducing kernel Hilbert space for n itself, its native 
space J\f K , and the second is the semi-inner product for the Hilbert space modulo Uj, M K j. It is 
important to note that norms for W™ and M K are equivalent, 



|w 2 m (M)- (4-5) 
This follows easily from kj ~ c^\~ m ~ c^K m j, for j large. 
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4.1 Interpolation with poly harmonic kernels 

If we choose kj > for j G J , the kernel will be positive definite, and so we can interpolate any 
/ G C(M). The form of the interpolant will be 

Ixf = J2 a Z< x >0- (4-6) 

with the coefficients being obtained by inverting f\x = Aa, where A = an d the 

components of o are the a^'s. 

Let's look at error estimates, given that / G iy™(M). Before we start, we begin with the 
observation that (|4.5l) and the usual minimization properties of || • || K imply that 



\\Ixf - fWwpM ^ c ^f ~ /IU < C||/|U < Call/Hw-CM)- 



This inequality was the key to proving Proposition |3.4[ where the reproducing kernel Hilbert space 
was itself W™(M). Hence, repeating the proof of that proposition, with the inequality changed 
appropriately, we obtain the estimate 

\\Ixf - /|U 2 (M) < Cm/II^M), / e W?(M), (4.7) 



which holds for a polyharmonic kernel (4.1), provided the conditions on X in Proposition 3.4 are 
satisfied. There are no restrictions on the two-point homogeneous manifold M. 

There are, however, restrictions on the smoothness of / - namely, / must be at least as smooth 
as k. Put another way, / has to be in the native space J\f K . 

In the case of § n or MP n , it is possible to remove this restriction and "escape" from the native 



space of k. We will first deal with From (4.2), we see that polyharmonic kernels on §>™ 
are spherical basis functions (SBFs), provided coefficients with j G J are taken to be positive. 
Moreover, kj ~ Xj m for j large. We may then apply [33j Theorem 5.5], with (ft, r, f3, [i — > <£, m, fj,, 0, 



to obtain the estimate below. (Note that for an integer k, the space Hk(S n ) in 33 is W§ (§ n ) here.) 
This result applies to functions not smooth enough enough to be in W™. 



Proposition 4.2 (Case ofS"). Let k be a polyhamonic kernel of the form (4-2), with deg(7r(a;)) = 
m, and with kj, j G J, chosen to be positive. Assume that the conditions on X C S n in Proposi- 
tion 3.4 are satisfied. If fi is an integer satisfying m > fj, > n/2 and if f G Wf(S n ), then, provided 
h is sufficiently small, 

\\Ixf-f\\ L2m <Ch»\\f\\ w » (§n) . (4.8) 

The n-dimensional, real projective space, MP™, is the sphere S n with antipodal points identified. 
Thus each x G MP™ corresponds to {x, —x} on S n . We will use this correspondence to lift the entire 
problem to S n , which is an idea used in |15[|21| in connection with approximation on SO(3). As we 
mentioned earlier, the distance d(x,y) = c?Kp™(x,y) has been normalized so that the diameter of 
MP n is 7T, rather than the ir/2 one would expect from its being regarded as a hemisphere of S n . In 
fact, the two distances are proportional to each other. The natural geodesic distance on projective 
space is simply the angle between two lines passing through {x, —x} and {y, —y}, with x, y G § n , 
which is just 8 = arccos(|x ■ y\), < < tt/2. The distance d|p», in which the diameter of M is it, 



satisfies dw n (x, y) = 2 arccos(|x • y|). If we use this in (4.2 ), then t = \x ■ y\ and k(x, y) = <I>(|x • y 



One can take this farther. In the case of MF n , the series representation for $ is 

*(t) = ^^4 2 • 2 \t). 

j=0 
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Since the Jacobi polynomials P n 2 ' 2 (t) are even or odd, depending on whether n is even or 
odd, and since the series for $ contains only polynomials with n even, it follows that $(— t) = <£(£). 
As a consequence we have that n(x,y) = • y). Thus, if we regard the variables x and y to be 
points on §™, the kernel k is even in both variables, and is nonnegative definite on §>™. 

Our aim now is to lift k to a polyharmonic kernel k* on §>™, with a view toward lifting the 
whole interpolation problem to §™. To begin, note that the eignvalues for — A§n have the form 

= j(j ' + n ~ 1) an d so those the MP™ may be written as Xj = AJjj. Next, let n(x) E II m be the 



polynomial associated with the kernel k in (4.1) and let k* = vr(A*) , except where tt(A*) < 0. 
For these we only require that the corresponding k* > 0. Define the function 

'^~ ) ( n ~2 n — 2 \ 

i=o 

which is associated with the polyharmonic kernel k*(x, y) = $>*(x- y) on §™. In constructing <!>* we 
have simply added an odd function to Thus, 

*(t) = !(<&*(*) + **(-*))• 

Start with the centers X on MP™. Each center in X corresponds to {£, — ^} on §>™, so we may 



lift X to X* = {±£ £ §>™ : {£, — £} 6 X}. Following the proof in 15, Theorem 1], we may show 
that qx' = \qx and hx' = hhx, where dw^ is used for q and h in X. In addition, we may, and 
will, identify each / 6 C(MP™) with an even function in C(§"). Now, interpolate an even / on §™ 
at the centers in X* . This gives us 

ix-m= £ ( 4 - 9 ) 

Note that k*(-x,C) = **((-a?) • C) = • (~C)) = «*O,-0- Since both C and -£ are in X*. 
it follows that Ix*f(~x) = £ Ce x* a*_ ( K*(x,0- Moreover, Ix*/(-£) = /(-£) = /(C) = 
Since interpolation is unique, and the two linear combinations of kernels agree on X* , they are 
equal, and so we have Ix*f{—x) = Ix*f{x) and a*_^ = a^. Because Ix*f is even, we have, 



from (4.9), that Ix*f(x) = Yy(ex* a (l{ K *( x ,C) + -())■ Since k*(x,() = <&*(x ■ (), we have 
K*(x, Q + K*(x, —0 = 2k(x, Q. Consequently, 

Ix*f(x) = a ( K (. x ,() = ^ac K ( x 'C), 
(ex* (ex 

where = a£ + a 4 ^ = 2a£. The sum on the right above is the interpolant to / on MP™, Ixf{x). 
Thus, I x *f(x) = I x f(x). 



Corollary 4.3 (Case of MP™). Let k be a polyhamonic kernel of the form (4-%)> with deg(7r(x)) 



m, and with kj, j £ J , chosen to be positive. Assume that the conditions on X C MP™ in 



Proposition 3.4 are satisfied. If \x is an integer satisfying m > fx > n/2 and if f £ W^MP 71 ), then, 



provided h is sufficiently small, 

\\Ixf - /||l 2( rp») < Ch»\\f\\ w » (wn) . (4.10) 

Proof. The metric dsKP n we are using for MP™ is twice the metric inherited from the sphere; that 
is, dsKP" = 2(fs§n . From the point of view of integration, it is easy to see that d/i^F" = 2 n / 2 d/i§« . 
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Thus various integrals and norms are merely changed by d-de pendent constant multiples - e.g., 
\\Ixf ~ /IU 2 (MP») = 2 1 - d / 2 ||/^/ - f\\ L2i &ny By Proposition 



4.2 



\!*xf - /lli 2 (s«) <Cih x *\\I x *f - fWw%(§ n ) < C 2 h x \\I x f - f\\ w g(wp n )- 



The inequality (4.10) then follows immediately from this and the remarks above. 



□ 



4.2 Error estimates for interpolation reproducing II j 

The interpolation operator associated with k that reproduces Hj = spanj^^ : j G J, 1 < s < dj} 



is 



Ix,jf = ^2 a^jK(x,i) +pj, pj G Uj, and ^ a^j(pj,s(0 = °> i G 

We wish to estimate the L2-norm of Ix,jf — /• To do that, we will need to obtain various equations 
relating Ix,jf, Ixf, and pj. (Since the choice of kj, j G J , obviously has no effect on Ix,j, we 
can and will assume that kj > 0, j G J .) Let Pj be the L 2 orthogonal projection onto Uj. Note 
that the interpolant Ix,jf = X^ex a £,J K i x i + PJ consists of two terms. It is easy to show that 
the first term belongs to (Tlj) 1 ~ and, of course, pj £ lij. Consequently, 



Pj(Ix,jf-f)=Pj-Pjf. 



(4.11) 



Next, because both interpolate /, the difference ^§ex( a £ ~~ a £,j) K ( x iO interpolates pj. This 
implies that 

Ixf - Ix,jf = IxPj - PJ, 



From this, we have 



\Ix,jf - /IU 2 (m) < \\Ixf - /IIl 2 (m) + Wpj ~ IxPj\\l 2 (m) 



(4.12) 



Since pj — IxPj\x = , the second term on the right can be estimated via Theorem 3.3, provided 
h is sufficiently small. Making the estimate yields this bound: 

Wpj ~ IxPj\\l 2 (m) < Ch m \\pj - IxPj\\w^(M)- 

In addition, by ( |4.5| ), the norm induced by k is equivalent to the VF^M) norm, and thus we have 
that \\pj — IxPjW w 2 m (M) < C\\pj — IxPjWk < CHpjIU; where the rightmost inequality follows from 
the minimization properties of IxPj in the norm || • From the definition of the k norm, where 
we have assumed kj = 1, j G J, it is easy to see that Hpj-Hk = ||pj r ||L 2 (M)- Combining the various 
inequalities then yields 

Wpj - Ixpj\\l 2 (M) < Ch m Wpj\\L 2 (M)- 

Using this on the right in ( 4.12[ ) gives us 



\Ix,jf ~ /||t 2 (M) < \\Ixf ~ /IIl 2 (m) + Ch m Wpj\\L 2 (M)- 



(4.13) 



Furthermore, employing (4.11) in conjunction with this inequality, we see that 

Wlx,jf ~ /|U 2 (M) < Uxf ~ /IU 2 (M) + Ch m {Wl X ,jf ~ /b 2 (M) + WPjfWleM) 

Choosing h so small that Ch m < | and then manipulating the expression above, we obtain this 
result. 
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Lemma 4.4. Let f G L2(M). For all h sufficiently small, we have 



\\Ix,jf - fh 2 (M) < 2 W T xf ~ fh 2 (M) + Ch m \\Pjf\\ L2[M) . (4.14) 
We now arrive at the error estimates we seek. 

Theorem 4.5. Let M be a two-point homogeneous manifold and let k be a polyhamonic kernel of 



the form (4-2), with deg(7r(x)) = m. Assume that the conditions on X C M in Proposition 3.4 



are 



satisfied. If fx is an integer satisfying m > p > n/2, then, provided h is sufficiently small, 



\Ix,jf - /IU 2 (M) < 




W 2 m (M)> 



H/ 2 M (M)' 



all M, and f G W^(M), 

M = § n ,MP n , and f G W£(M). 



(4.15) 



Proof. In (4.14), we have \\Pjf\\ L2 (M) < \\fh 2 (U) < 



\\Ixf ~ /IU 2 (M) follow from (4.7), (4.8), and (4.10). Combining these yields (4.15). 



VK 2 M (M) 



; in addition, the various bounds on 

□ 



4.2.1 Optimal convergence rates for interpolation via surface splines 



For spheres, the interpolants constructed from the restricted thin-plate splines defined in (4.3) 
converge at double the rates discussed above - namely, 0(h 2m ) rather than 0(h m ) -, provided / 



is smooth enough. This also applies in SO(3) for the kernels given in (4.4). The spaces Hj being 



reproduced here are, for E> n , the spherical harmonics of degree m — 1 or less; and for 50(3), the 
Wigner D-functions for representations of the rotation group, again having order m — 1 or less. The 
precise result is stated below. 



Proposition 4.6 ( 20, Corollaries 5.9 & 5.10]). Suppose m > n/2 and let Ix,j denote the inter- 
polation operator corresponding to the restricted surface spline defined in (4-3) for s = m — n/2. 



Then, there is a constant C depending on p, m and n such that, for a sufficiently dense set X C S n ; 
and for f G C 2m (§ n ), the following holds: 



\Ix,jf ~f \\oo <Ch 



2/ii 



C 2m (S")> 



Similarly, for the kernels in (44), if f G C 2m (SO(3)), then \\Ix,jf - f\\oo < Ch 



-,2m 



C 2m (SO(3))- 



4.3 Lagrange functions and Lebesgue constants 



If x%( x ) is the Lagrange function associated with a kernel k and a space Hj, then, by (2.7), the 
weights in the quadrature formula have the form eg = j M x^(x)dp(x), £ G X. Our aim is to use 
properties of x§ to obtain bounds on these weights. Before we do this, we will need the following 
decay estimates for xs,- 



Proposition 4.7 ( |20, Theorem 5.5]). Suppose that M is an n-dimensional compact, two-point 
homogeneous manifold and that k is a polyharmonic kernel, with deg(7r(»)) = m, where m > n/2. 
There exist positive constants ho, v and C, depending only on m, M and the operator C m = ir(— A) 
so that if the set of centers X is quasi-uniform with mesh ratio p and has mesh norm h < ho, then 
the Lagrange functions for interpolation by k with auxiliary space Uj satisfy 



\xd*)\<c P m - n/2 



max I exp 



d(x,m ,h 



,2m 



(4.16) 
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There are two results that we want to obtain. First, we want to estimate the size of each weight 
in terms of the mesh norm h; and second, we want to estimate X^gex \ c %\- 

We only need to estimate ||x^IUi(M)> since |cg| < ||xclUi(M)- Ovrr approach is to divide the 
manifold into a ball B = B^R h , having radius Rh and center £, and its complement B^. The radius 
Rh is the "break-even" distance in (4.16); it is obtained by solving exp (— j^Rh) = h 2m for Rh. The 
result is R h = ^fh\ \og{h)\. By (|4.16|)7we obtain 



\xt(x)\dfi(x) < Cp m - n/2 fi(M)h 



2 m 



Next, again by (4.16), we have that 



Xt(x)\d»(x) < C'p m - n / 2 

B JO 



-™/h r n-l dr < C > p m-n/2 



-rv/h r n-l dr 



(n-l)\{h/v) r - 



Consequently, |cg| < ||xclUi(M) < Cp m n l 2 {hjv) n {\ + h 2m n v n ). Because 2m > n, for h small 
enough, it follows that 

k| < ||xelk(M) < C'p m - n ' 2 h n , (4.17) 

where C = C"(M, k, J). 

The next result concerns the boundedness of the Lebesgue constant, which played a significant 
role in section |2.3| in the analysis of the stability of the quadrature operator. 



Proposition 4.8 ( [20| Theorem 5.6]). Let the notation be the same as that in Proposition ^.7. If 
h < ho, then the Lebesgue constant, Ax t j, K = max xeM Sgex 1x^(^)1? associated with X, J , and n, 
satisfies the bound Ax,j,k < C, where C depends only on m, p, and M. 



4.4 Quadrature via polyharmonic kernels 

In this section, we will employ the various properties of polyharmonic kernels, which are of course 
invariant, to obtain results concerning accuracy and stability of the associated quadrature formulas. 

Theorem 4.9. Suppose that M is an n- dimensional compact, two-point homogeneous manifold 
and that k is a polyharmonic kernel, with deg(7r(x)) = m, where m > re/2. Let X C M be 
a finite set having mesh ratio px < p, Hj = span{0 JjS : j : £ J , 1 < s < dj}, and Qv x be 



the corresponding quadrature operator given in Definition 2.3 . The norm of Qv x is bounded, 
WQvx llc(M) — MM)Ax,j\k — C mjPj M, and the error satisfies the estimates 



\QvAf) 



fdfil < 




wr(M) , all M, and f G W 2 m (M), 

w g(M), M = S",IRP n , and f £ W£(M), n/2 < p < m. 



(4.18) 



Finally, the standard deviation o~q from Proposition 2.4 satisfies o~q < Ca v h n l 2 , where C 
C(p,m,M,J). 



Proof. The norm estimate follows from (2.16) and Proposition 4.8 The error estimate is a conse- 



quence of (2.14) and Theorem 4.5 Finally, the bound on <tq is a consequence of Proposition 
Proposition 4.8, and the bound on ||x^|| in (4.17). 



2.4 



□ 
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At present, the most important compact two-point homogeneous manifolds are spheres and pro- 
jective spaces (especially S 2 and SO (3)). As we discussed in section 4.2.1 the restricted thin-plate 



splines (4.3) on § n and similar kernels (4.4) on SO(3) give interpolants with optimal convergence 



for smooth target functions. This also is reflected in the accuracy of the corresponding quadrature 
formulas: 



Corollary 4.10. Let X, p x , p, m,H.j be as in Theorem]^ Take M = S n or SO(3) and Q to 
be the quadrature operator corresponding to the restricted surface spline defined in |^.3| ) or to that 
for the polyharmonic kernel in (4-4)- If X is a sufficiently dense in M, then there is a constant 
C = C(m, n, p) such that for a sufficiently dense set Ic§" and for all f 6 C 2m (S n ), 

\Q(f)~ ! fdp\<Ch 2m \\f\\ c2m{sn) , where C = C(m,n,p), 
Jm 



in addition to the bounds on \\Q\\ and o~q from Theorem 4-9 holding. 



Because of applications in physical sciences and engineering, § 2 is undoubtedly the most im- 
portant of the manifolds treated here. In the next section, we give some numerical examples for S 2 
and the m = 2 surface spline <£(t) = (1 — t) log(l — t) (i.e. the thin plate spline r 2 log(r 2 ) restricted 
to §> 2 ) that validate the above theory. 



5 Numerical results for S 2 

We begin with a brief overview of how the quadrature weights can be computed in an efficient 
manner for the m = 2 surface spline <J>(i) = (1 — t) log(l — t) using the local Lagrange preconditioner 
developed in [TlJ. This is followed by a description of the nodes used in the numerical experiments 
and some properties of the resulting surface spline quadrature weights and their stability. Finally, 
we give some results validating the error estimates from the previous section. 



5.1 Computing the quadrature weights 



The m = 2 restricted surface spline kernel is conditionally positive definite of order 1 and the 
finite dimensional subspace II associated with it consists of all spherical harmonics of degree < 1, 
i.e. IT := span{Yo,0) Yl,o, ^1,1, ^1,2}) where Y^fc is the degree I and order k spherical harmonic and 
< k < 21 + 1. Given a set X = {xj}^ =l of distinct nodes on S 2 , the quadrature weights c for 
this kernel can be computed by first solving (2.10) for cj_ and then computing c via (2.11). In 
these equations, Ajj = (1 — Xi ■ Xj) log(l — xi ■ Xj), i,j = 1, . . . , N, and ^ is the A^-by-4 matrix 
with columns = io,o (zi) and ^i t k+2 = ^i,fe(aJi), for i = 1, • • • , N, k = 0, 1, 2. Additionally, 

J = [4vr 0] T and J = 2vr(41og(2) - 1). 

Since ^ only has 4 columns, computing ^(^ T ^ / )~ 1 J in (2.10) can be done rapidly using, 
for example, QR decomposition. Thus, the bulk of the computational effort in computing the 
quadrature weights is in solving for c± in (2.10). Since the matrix A is dense, direct methods cannot 
be realistically applied for large N, and one must then resort to iterative methods. However, for 
iterative methods to be useful, one must apply an effective preconditioner to the system. In [11] , 
we developed a powerful preconditioner for ( 2.10| ) based on local Lagrange functions and combined 
it with the generalized minimum residual (GMRES) iterative method 36 . The basic idea of the 
preconditioner is, for every node xj £ X, to compute the surface spline interpolation weights for a 
small subset of nodes about Xj consisting its p = M(logN) 2 nearest neighbors, where M is suitably 
chosen constant. The data for each interpolant is taken to be cardinal about Xj. This is similar to 
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the preconditioner used in [8] for interpolation on a 2-D plane, however, in that study the number 
of nodes in the local interpolants did not grow with N and it was observed that the preconditioner 
broke down as N increased. As demonstrated in by allowing the nodes to grow very slowly with 
N the preconditioner remained effective and the total number of iterations required by GMRES to 
reach a desired tolerance did not increase with N. We refer the reader to for complete details 
on the construction of the preconditioner. 

In the numerical results that follow, we used the preconditioned GMRES technique from 
with p = 2[~(log iV) 2 ] and a relative tolerance of 10 -12 to solve for c± in (2.10). We then used this 
in (2.11) to find c. Table [2] lists the number of GMRES iterations required to compute c± for three 
different quasi-uniform node families, which are described in the next section. We see from the 
table that the number of iterations remains fairly for constant as N grows for all three families of 
nodes. 



Icosahedral 


Fib 


onacci 


Min. 


Energy 


N 


Iterations 


N 


Iterations 


N 


Iterations 


2562 


8 


2501 


9 


2500 


9 


10242 


7 


10001 


8 


10000 


8 


23042 


7 


22501 


11 


22500 


7 


40962 


7 


40001 


8 


40000 


8 


92162 


6 


62501 


10 


62500 


7 


163842 


7 


90001 


10 


90000 


7 


256002 


6 


160001 


8 


160000 


7 


655362 


6 


250001 


10 


250000 


7 



Table 2: Number of GMRES iterations required to compute c± in (2.10) using the preconditioned 
iterative method developed 111] for determining the m = 2 surface spline. The quadrature weights 



are then computed from (2.11). In all cases, the relative tolerance of the GMRES method was set 
to 10- 12 . 



We conclude by noting that as part of the iterative method, matrix-vector products involving 
A must be computed. Since A is dense, this requires 0(N 2 ) operations per matrix-vector product. 
In the computations performed for this study, these products were computed directly, making the 
overall cost of the weight computation 0(N 2 ). In a follow up study we will explore the use of fast, 
approximate matrix- vector products using the algorithm described in [28] . By using this algorithm, 
it may be possible to reduce the total cost of computing the weights (or a surface spline interpolant) 
to O(NlogN). 



5.2 Nodes, weights, and stability 

We consider three quasi-uniform families of nodes for the numerical experiments. The first is the 
icosahedral nodes, which are obtained from successive refinement of the 20 spherical triangular faces 
formed from the icosahedron. The second are the Fibonacci (or phyllotaxis) nodes, which mimic 
certain plant behavior in nature (see, for example, |14| and the references therein). The third are 
the quasi-minimum energy nodes, which are obtained by arranging the nodes so that their Riesz 
energy is near minimal |23| . In the examples below, a power of 3 was used in the Riesz energy and 
the nodes were generated using the technique described in [5] . The mesh norm for all three of these 
families satisfies h ~ -7=, where N is the total number of nodes. The mesh ratio p stays roughly 
constant for the Fibonacci and quasi-minimum energy nodes as N increases. For the icosahedral 



21 




x 10" 
I 6.5 

-6 

5.5 

-5 

-4.5 

-4 

I 3.5 



x 10" 

5.9 
5.8 
5.7 
5.6 
5.5 
5.4 
5.3 



(a) Icosahedral nodes and weights, N = 23042 (b) Fibonacci nodes and weights, N = 22501 



x 10" 
6.2 



5.8 



5.6 



5.4 



5.2 



(c) Quasi-minimum energy nodes and weights, N = 22500 

Figure 1: Visualization of the different families of nodes used in the numerical experiments and the 
corresponding quadrature weights for the m = 2 spherical spline. The nodes are plotted on § 2 using 
an orthographic projection about the north pole. Each node has been given a color corresponding 
to the value of the quadrature weight for that node. 



points p grows slowly with N since the spacing of the nodes decreases faster towards the vertices of 
the triangles of the base icosahedron than at the centers [37]. However, in the numerical examples 
that follow, this increase seems to be of little concern. All three of these families of nodes are quite 
popular in applications; see, for example [13 29 35,41 for the icosahedral nodes, (27 39 42 for the 



Fibonacci nodes, and [9j[ToJ[38j[48] for the quasi-minimum energy nodes. Quadrature over these 
nodes also plays an important role in applications, for example, for computing the mass of a certain 
quantity, the energy for a certain process, or the spectral decomposition of some data. 

It should be noted that previous studies have been devoted to developing quadrature formulas 
and error estimates for icosahedral pi Ch. 5] and Fibonacci [22] nodes. However, these results rely 
on the specific construction of the node sets and cannot be applied to more general quasi-uniformly 
distributed nodes such as the quasi- minimum energy nodes. 
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Figures [T] (a)-(c) show examples of the different families of nodes and provide a visualization 
of the values of the corresponding quadrature weights for the m = 2 surface spline. The geometric 
pattern of the icosahedral nodes in part (a) of this figure are clearly reflected in the values of the 
corresponding quadrature weights. This is also true of the Fibonacci nodes in part (b), which have 
a slight clustering near their seed value (in this case the north pole), but then are quasi- uniformly 
distributed. There are no clear patterns for the minimum energy nodes in part (c) of Figure [TJ as 
these are not distributed in a discernible pattern. Looking at the color bars in each of the plots 
we see that the range of values of the weights is similar between the different node families and 
comparable to 1/N, and are also positive. 



Mill. Energy 

IcoHjihedral 

Fibonacci 

-jy-1/2 



Figure 2: Estimate of the standard deviation oq in (2.17) for the different families of nodes. The 
estimate is based on sample sets of 500 iV-point quadratures of different independent, identically, 
distributed, zero mean (quasi) random data with a standard deviation of 1. This confirms the 



stability estimate in the last part of Theorem 4.9 



To estimate oq in (2.17), and hence the stability of the quadrature weights in the presence of 
noise, we performed the following experiments. For each family of nodes, a value of N was selected 
and then a sample set was generated consisting of 500 values of the iV-point quadrature of different 
independent, identically, distributed, zero mean (quasi) random data. The standard deviation of 
the sample sets was then computed to estimate <jq. The results are plotted in Figure [2] as a function 
of ./V for each of the three families of nodes. Comparing the results to the dashed line on the plot, 
we see that the estimated ctq decreases like 0{N~ l l 2 ), or like 0(h) since h ~ N~ 1 / 2 for these 
families of nodes. This is in perfect agreement with the rate predicted by the last part of Theorem 



4.9 for the surface splines. 



5.3 Convergence results 



Two target functions of different smoothness were used to test the error estimates of Theorem 4.9 
and Corollary 4.10 The target functions were chosen so that the Funk-Hecke formula (see, for 
example, j5J §2.5]) could be used to determine their exact integral over S 2 . Letting x,x c £ E 2 and 
g be a zonal kernel, i.e. g(x, x c i) = g(x ■ x c ), such that g G L [— 1, 1], the Funk-Hecke formula gives 
the following result: 



/ 



g(x ■ x c )Yg i k(x)dfJ l (x) 



\-KCLl 
21 + 1 



Yi,k( x c 



(5.1) 
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where Y^f- is any degree £, order k spherical harmonic, and a? is the £th coefficient in the Legendre 
expansion of g(t). 

The following two kernels were used in constructing the target functions: 

gi{t ) = -(2-2t) l /\ (5.2) 

92(t) = - 1 ~ £ -3 (0<e<l), (5.3) 

(1 + e 2 -2et)2 

which are known as the potential spline kernel of order 1/4 and the Poisson kernel, respectively. 
The Legendre expansion coefficients for these kernels are given as follows (see (4l): 

(-l)Wy/2(r(f)) 2 (2l+l) 

91 '■ ai= r(l-i)ra + i) ■ (5 - 4) 

ff2 : a, = {21 + l)e«. (5.5) 

The smoothness of these kernels is of course determined by the decay rate of the Legendre coeffi- 
cients ag. For <7i, we have d£ ~ £~ 3 / 2 , which means g\ belongs to every Sobolev space W^iS 2 ) with 
/j, < §. While for g 2 the Legendre coefficients decay exponentially fast, which means 52 £ C 00 ^ 2 ), 
analytic in fact. 

Using the above results, we define the following two target integrands: 

41 



A 0*0 = ^ sign(y 2 o ifc (x c ))y 2 o )fc (x)ffi(% • x c ), (5.6) 
fc=i 

41 

M^) = ^ sign(l2o,fe(^c))i20,fe(a;)52(^ ■ x c ), (5.7) 



fc=i 



where x c = (cos(-2. 0281) sin(0.76102), sin(-2. 0281) cos(0.76102), sin(0. 76102)) and e = 2/3 for g 2 ; 
see Figure [3] (a) and (b) for plots of these respective functions. Integrating these functions over S 2 



24 



and applying the Funk-Hecke formula (5.1) gives 



fi(x)dn(x) = 0.014830900415995, 
f 2 (x)dfi{x) = 0.032409262543520. 



(5.8) 
(5.9) 



Both /i and f 2 inherit their smoothness directly from g\ and g 2 , respectively. Thus, f\ belongs to 
every Sobolev space W^(S 2 ) with fi < § and f 2 G C°°{§ 2 ). 

Tables [3] and [4] display the relative errors in the iV-point quadrature of fi and f 2 , respectively, for 
the different family of nodes, while Figures |4]^a) and (b) display these respective results graphically 
on a log — log scale. Focusing first on the results for f± in Figure [4]^a) and comparing the results to 
the included dashed line, it is clear that for all three families the error is decreasing approximately 
like 0(N~ 1 ' 25 ). Since h ~ iV" 1 / 2 for these families of nodes, this observed rate of decrease in the 
error is approximately 0(/i 2,5 ), which is precisely the rate predicted by the estimates in Theorem 
4.9| for functions with fi's smoothness. Doing a similar comparison of the results for f 2 in Figure 
4[b), we see that the errors associated with the icosahedral nodes very clearly decrease like 0(N~ 2 ). 
The results for the other two nodes are not as clear, but for larger values of N the errors do appear 
to be decreasing approximately like 0(N~ 2 ). Again, because of the relationship between h and N 
for these nodes, the errors for f 2 are thus decreasing approximately like 0(h ). This is the expected 
rate from Corollary 4.10 since f 2 is infinitely smooth and we have used the m = 2 surface spline. 
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N 


Rel. Error 


N 


Rel. Error 


N 


Rel. Error 
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l 
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3 
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3 
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2 
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3 
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2 
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3 
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1.040 x 10" 


3 
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2 


92162 


2.273 x 10" 


3 


62501 


6.460 x 10" 


4 


62500 


6.989 x 10- 


3 


163842 


1.107 x 10- 


3 


90001 


4.068 x 10" 


-4 


90000 


4.393 x 10- 


3 


256002 


6.330 x 10- 


4 


160001 


1.844 x 10" 


-4 


160000 


2.083 x 10" 


3 


655362 


1.957 x 10- 


4 


250001 


1.085 x 10" 


-4 


250000 


1.228 x 10" 


3 



Table 3: Relative error in the iV-point quadrature of the "rough" target function f\ in (5.6) for the 
different families of nodes. 

We conclude by noting that the nodes and quadrature weights used in the numerical experiments 
above are available for download from |47| . 

6 Quadrature on Manifolds Diffeomorphic to Homogeneous Spaces 

Numerical integration of functions defined on smooth surfaces that are diffeomorphic to S 2 arise 
in a number of applications. For example, the shape of the earth, as well as the other planets, is 
a "flattened sphere" - i.e., an oblate spheroid, which is diffeomorphic to a sphere. To numerically 
compute integrals over the earth's surface then requires quadrature formulas over oblate spheroids. 
In addition, the need for numerical integration over various surfaces also arises in boundary element 
formulations of continuum problems in 1R 3 |2l Ch. 6]. In this section, we discuss quadrature in a 
more general context. We will show how to use quadrature weights for a homogeneous manifold 
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Table 4: Relative error in the iV-point quadrature of the "smooth" target function /2 in (5.7) for 
the different families of nodes. 



M to obtain an invariant, coordinate independent quadrature formula for a smooth Riemannian 
manifold L diffeomorphic to M. Of course, the case in which L is an oblate spheroid and M = S 2 
is of special interest. 

Recall that a C°° manifold L is diffeomorphic to M if there is a C°° bijection F : M — > L. Let 
gij be the Riemannian metric for M. Suppose that L has the metric hij. Consider a local chart 
(U, (f) : U — )■ W 1 ) near po £ M. Using this chart we have coordinates x = (x , . . . , x n ) = 4>{p) and a 
parametrization p = 

We can use (U, 4>) to produce a local chart (V, ip : V — > W 1 ) near qo = F(po) G L. Simply let 
V = F(U) C L and tp = <po F . The coordinates for V are thus x = tp(q), and so the two metrics 
gij and hij can be expressed in terms of the same set of coordinates. In these coordinates, the 
volume elements are given by 



dfiM(x) = &et(gij(x)) dx 1 ■ ■ ■ dx n and dfij^(x) = det(/ijj(x) dx 1 ■ ■ ■ dx n . 
It follows that 

dfJ,h(x) -- 



' det(hij(x)) 



det(gij(x)) 
j 

w(x) 



dnm(x). 



(6.1) 



Suppose that we now make a change of coordinates, from x = <p(p) to new coordinates y = tp(p), 
or x = x(y). Let x'{y) be the Jacobian matrix for the transformation, and let J(y) = det(x'(y)). 
In y coordinates, the metrics are g{y) = (x' (y)) T g(x(y))x' (y) and h(y) = (x'(y)) T h(x(y))x'(y). 
Consequently, we have that 

det{gij(y)) = J{y) 2 det(gy («(?/))) and det(hij(y)) = J(y) 2 det (hi j(x{y))), 



and, furthermore, that 



w{x{y)) 



l det{hjj{x{y))) _ / dggggj 
det( gij (x(y))) V det(gUy)) 



w(y). 



This means that w o <j)(p) = w o ip(p) =: W(p) is thus a scalar invariant that is independent of the 
choice of coordinates. In terms of integrals, we have 



foF(p)W(p)d m (p). 



(6.2) 
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(a) Rough target, f\. 



(b) Smooth target, f2- 



Figure 4: Relative errors in the iV-point quadratures of (a) f\ and (b) f2 for the three different 
families of nodes. The mesh-norm for all three families of nodes satisfies h ~ iV -1 / 2 , so that the 
dashed lines in the figures indicate convergence like (a) 0(h 2 ^) and (b) 0(h 4 ). These are the 
predicted convergence rates from Theorem 4.9 and Corollary |4. 10 respectively. 



An invariant, coordinate independent quadrature formula for L can be obtained from the one 
for M. We have the following result. 

Proposition 6.1. Let X denote the set of centers on L and let X' = F~ l (X) be the corresponding 
set on M. Suppose that we have a quadrature formula for M with weights {C^}^ e x' ■ Then we 
have the following quadrature formula for L, 

QUf) ■= E where c « := wiF-^o^F-no- ( 6 - 3 ) 

Proof. Applying the quadrature formula for the homogeneous manifold M to the integral on the 



right-hand side of (6.2) yields 



Qm{f°F(p)W(p))= /°*X£'WOQ'- 
Since £' = .F _1 (£), we have /oF(0 = /(£) and W{£') = W{F~ 1 {^)). Taking c 5 = W(F~ 1 {£))C F -i 



(0 



we obtain (6.3). □ 



Oblate spheroid Consider the oblate spheroid L, x 2 + y 2 + z 2 /a 2 = 1, where < a < 1, and 
the 2-sphere S 2 (= M), X 2 + Y 2 + Z 2 = 1. The diffeormorphism between the two manifolds is 
F(X,Y,Z) := (X,Y,aZ); that is, (x,y,z) = (X,Y,aZ). Our aim is to find the scale factor W. 
Since the end result will be coordinate independent, we choose to work in spherical coordinates 
(6, 4>), where 8 is the longitude and <j) is the latitude on S 2 . (The north pole is (0, 0, 1).) Obviously, 
for L we have x = sin# cos0, y = sin 9 sin <f>, z = a cos 6. The metric for the sphere is dS 2 = 
d6 2 + sin (9)d(f) 2 . The metric for L is the Euclidean metric dx 2 + dy 2 + dz 2 on M 3 restricted to L. 
Making a straightforward calculation, one can show that the metric for L is 



ds 2 = (cos 2 6 + a 2 sin 2 9)d9 2 + sin 2 9d(j) 2 = (a 2 + (1 - a 2 ) cos 2 9)d9 2 + sin 2 
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Consequently, the volume element for L is 

(i/iL = \/ a 2 + (1 — a 2 ) cos 2 9 sin 6 d6d(j) = \J a 2 + (1 — a 2 ) cos 2 9 d//§2 

v v ' 

To put this in invariant form on §> 2 , we use Z = cos 9 to obtain PF(A, Y, Z) = a 2 + (1 — a 2 )Z 2 . 
Pulling this back to L, we have W o F~ 1 (x,y, z) = \Ja 2 + (a~ 2 — l)z 2 . The weights for Ql are 
thus 

c ? = y/a? + (a- 2 -l)g 

As we mentioned above, the earth is approximately an oblate spheroid. The parameter a is the 
ratio of the polar radius to the equatorial radius. The flattening of the earth is / = 1 — a, and has 
the approximate value / ~ 1/300 (cf. Earth Fact Sheet |46|). From this, we get that a ~ 299/300, 
and so W o F~ 1 (x, y, z) ~ V0.993 + 0.007.2 2 . This factor varies between 0.9967 and 1, about 0.3%. 
Even so, it could affect the accuracy of the quadrature formula for functions with large values near 
the equator. As an aside, Jupiter has / « 1/15 (cf. Jupiter Fact Sheet (ellipticity) [46]), and for it 
the change in W o F~ l would be a hefty 7%. 
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